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Abstract 

Although Nose’s thermostated mechanics is formally consistent with Gibbs’ canonical ensem¬ 
ble, the thermostated Nose-Hoover ( harmonic ) oscillator, with its mean kinetic temperature 
controlled, is far from ergodic. Much of its phase space is occupied by regular conservative tori. 
Oscillator ergodicity has previously been achieved by controlling two oscillator moments with 
two thermostat variables. Here we use computerized searches in conjunction with visualization 
to find singZy-thermostated motion equations for the oscillator which are consistent with Gibbs’ 
canonical distribution. Such models are the simplest able to bridge the gap between Gibbs’ 
statistical ensembles and Newtonian single-particle dynamics. 
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I. ERGODICITY OF THE EQUATIONS OF MOTION 


Gibbs’ statistical mechanics is based on summing contributions from ensembles of sim¬ 
ilar systems. His microcanonical ensemble includes all the states of a given system which 
have the same energy. These energy states are accessible to a single “ergodic” system 
obeying Newtonian mechanics^. A periodic hard-disk or hard-sphere fluid is the usual 
example. Comparisons of Monte Carlo microcanonical-ensemble averages with molecular 
dynamics dynamical averages have conhrmed this equivalence, even for small systems of 
just a few particles^. 

Certainly Boltzmann and Gibbs both realized that all states need to be accessible to the 
dynamics in order for the dynamical and phase averages to correspond. The Ehrenfests 
had a practical dehnition of “quasiergodicity”. They used the word to indicate that the 
dynamics eventually comes “arbitrarily close” to all states. Their idea expresses very well 
our own view of what we call “ergodicity” in the present work. 

Gibbs’ canonical ensemble sums Boltzmann-weighted contributions from all energy 
states. The underlying idea is that the system of interest is weakly coupled to a heat 
reservoir with an ideal-gas density of states characteristic of a hxed kinetic temperature 
T . Shuichi Nose^i^ developed a dynamics consistent with the canonical distribution by 
including a “time-scaling” variable s and its conjugate momentum ps in the equations 
of motion. The new momentum ps acts as a thermostat variable capable of exchanging 
energy between the system and a heat reservoir at temperature T. 

Hoover showed that a harmonic oscillator thermostated in this way is not at all ergodio^. 
That is, there is no initial condition from which the oscillator is able to access all of 
its phase-space states. Instead, this thermostated oscillator has a nonergodic highly- 
complicated multi-part phase-space structure^. There are inhnitely-many regular non- 
chaotic orbits embedded in a single chaotic sea. Where “Chaos” controls the motion two 
closeby points, ri(t) and r 2 (t) , tend to separate from one another exponentially fast, 
either forward or backward in time. Such a motion is said to be “Lyapunov unstable”. 
The averaged separation rate is described by the largest Lyapunov exponent, Ai : 

6 = \r2-ri I ~ ; Ai = ( Ai(f) ) . 

The time-averaged Lyapunov exponent Ai is computed as an average of the instan¬ 
taneous local Lyapunov exponent, Xi{t). The local value is only rarely zero, even on 
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conservative tori, where the long-time averages vanish. We illustrate Ai(t) for the Nose- 
Hoover oscillator in Figure 1. We choose the simplest equations of motion, 

{ q = p ] P = -q - Cp \ C = p^-i}[nh]. 

They are time-reversible: any time-ordered sequence { +q, } satisfying the motion 

equations has a time-reversed backward twin, { -|-g, —p, —C, } satisfying the same equa¬ 
tions. The Nose-Hoover oscillator in { g, p, ^ } space is an improved and simplihed version 
of Nose’s dynamics, which occupies a four-dimensional { q,p,s,ps } space^i^. 

For this Nose-Hoover oscillator we have computed the local Lyapunov exponent on a 
grid of about a million points by the simple expedient of integrating backward in time and 
then forward, for a time of 100 in both directions. The “reversed” trajectory going forward 
in time can be compared to a nearby constrained “satellite” trajectory. We compute the 
instantaneous value of the time-dependent Lyapunov exponent just as the C = 0 plane 
is crossed. In the Figure red corresponds to the most positive exponent value and blue 
to the most negative. Within the chaotic sea the largest ( time-averaged ) Lyapunov 
exponent is O.OlSg . See Reference 6 for details. 

Outside the chaotic sea lie an inhnite number of regular orbits. All have a largest 
time-averaged Lyapunov exponent ( and also a smallest ) of zero. Because the oscillator 
is prototypical of systems with smooth minima in their energy surfaces a considerable 
effort has been made to hnd motion equations providing Gibbs’ canonical distribution for 

II. FEEDBACK CONTROL OF OSCILLATOR MOMENTS 

For simplicity we choose units of force, mass, time, and temperature corresponding to 
choosing the oscillator force constant, mass, angular velocity, and Boltzmann’s constant 
all equal to unity. In these units and without any thermostating the oscillator motion 
equation is q = p = —q . Because distribution functions for the oscillator’s displacement 
and momentum can be described in terms of their moments ( ^ ^ j^a^ural to 

control oscillator force and velocity moments with feedback variables such as ( and : 

{q = p-^q-, p = -q-Cp} ■ 
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Figure 1: This Nose-Hoover oscillator phase-space section corresponds to the plane C = 0 . The 
coloring reflects the local value of the instantaneous Lyapunov exponent at each { q,p,0 ) point, 
with red least stable and blue most. The distributions of{q,pX} in the chaotic sea are compared 
to Gibbs’ Gaussian distributions at the right. The white space indicates nonchaotic regions filled 
with regular nested tori, some of which are shown. In the chaotic sea Ai = ( Ai(t) ) = 0.0139 . 

The time dependence of the friction coefficients ({t) and ^{t) can be arranged so as to 
control one or more of the oscillator moments : 

C = (pVT)-1^(p^)=T; e = (gVT^) - 3(gVT) ^ ) = 3T( ) .... 

Bulgac and Kusnezov, along with their coworkers Bauer and considered a variety 

of simple systems and concluded that cubic contributions to the control equations, such 
as those in the [ HH ] and [ JB ] equations below, were especially useful in promoting 
chaos and ergodicity. 

Over thirty years dozens of investigators explored the ergodicity of thermostated 
oscillators^^—. Three successful models, the Hoover-Holiaiil^, Ju-Bulgaoi^, and Martyna- 
Klein-Tuckerman modelsl^ resulted. With all of the thermostat relaxation times set equal 
to unity, these models have the following forms : 

q = p; p=-q-Cp- {^pVT) ; C = (pVT) - 1 ; ^ = (pVT^) - ^{pVT) ; [ HH ] 
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q = p-, p=-q-C^p-{^ pVT) ; C = (pVT) - 1 ; ^ = (pVT^) - ^P^T) ; [ JB ] 
q = p; P = -q-Cp; C = {pVt)-i-^C; e = C'-i. [mkt] 

In all three cases the phase-space continuity equation, 

(df/dt) = -{dfq/dq) - (dfp/dp) - (dfC/dC) - {dfi/dO • 

shows that the motion equations are consistent with Gibbs’ canonical distribution for the 
{q,p) coordinate-momentum pair : 

The distributions of the control variables ( and ^ are likewise Gaussians : 

fuH = Imkt oc ; fjB oc . 

If the dynamics is ergodic, hlling out the full distributions, these models reproduce Gibbs’ 
canonical distribution. 

Although there are no theoretical proofs that these dynamics obey the phase-space 
distributions and relationships there is by now abundant numerical evidence that the 
three two-thermostat approaches given above are ergodio^^. Recent work, particularly 
that of Patra and Bhattacharyai^^i^i^ led us to search for even simpler models, with three 
equations rather than four, for thermostated oscillators. We describe our own quite novel 
and successful hndings next. 

III. SINGLY-THERMOSTATED ERGODIC OSCILLATOR MODELS 

A hrst look at the possibility of thermostating an oscillator with a single friction 
coefficient corresponds to a variety of separate models. We consider two of them here. 
Both include cubic functions as suggested by Bulgac, Kusnezov, Ju, and Bauer. The hrst 
oscillator model controls the huctuation of the kinetic energy : 

q = +p-p=-q- aiCY/T) ; ( = a[ {j? jTf - Z{j? jT) ] . 


The second controls the huctuation of the force : 

g = +P - /3(3g ; p = _g ; ( = /5[ [q^ jT) - 1 ] . 
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Figure 2: Single-thermostat cubic control enforcing the fourth-moment condition, { ) = 

3( ) with a and T both equal to unity and /3 = 0 . These choices show large gaps in the cross 

section where the time-averaged Lyapunov exponent vanishes. The probability densities within 
the chaotic sea are shown at the right. The initial condition used here and in all succeeding 
figures to access the chaotic sea is ( C ) = ( 0, 5,0 ) . The three horizontal nullclines at 
{ p = — \/3j 0, -|-\/3 } reflect the vanishing of the phase-space velocity component normal to the 
C = 0 plane. Ai = 0.1108 in the chaotic sea. 

Neither of these approaches is successful. In Figures 2 and 3 we show the chaotic seas 
corresponding to these two models with hrst a and then (3 equal to unity. In both cases 
we choose unit temperature, T = 1 . We plot {q-iP) points whenever the friction coefficient 
Q changes sign. These models both contain holes in the sea hlled with regular toroidal 
regions. The one-dimensional distribution functions shown at the right of these hgures 
give an alternative view of the models’ inability to describe Gibbs’ canonical distribution. 
Both of these single-thermostat models are failures. In addition to plotting cross sections 
and probability densities one can evaluate the likelihood of deviations from Gibbs’ values 
as measured by the statistic described in Wikipedia, Numerical Recipes, and many 
other texts. 
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Figure 3: Single-thermostat cubic control enforcing the second moment condition { ) = 1 

with /3 and T set equal to unity and a = 0. These choices show the presence of at least 
twenty gaps ( or “holes” ) in the cross section where the friction coefficient C vanishes. The 
one-dimensional probability densities within the chaotic sea are shown at the right along with 
the Gibbs’ distributions from the canonical ensemble. Ai = 0.0905 in the chaotic sea. 

At least in retrospect it is natural to consider the possibility that a single friction 
coefficient ( might somehow control two moments simultaneously, rather than just one. 
For example, consider simultaneous control of fluctuations in both the force ~ ( ) and 

the kinetic energy ~ ( — 3p^T ) : 

q = +p- ] P = -q- a(CV/T) ; 

c = f3[ (q^T) - 1 ] -f «[ (pVT)2 - 3(pVT) ] [ HS ] . 

Our numerical work indicates that this [ HS ] ( Hoover-Sprott ) idea has merit. With 
two free parameters { a, (3 ) there are an inhnite number of models which could be tested 
against the predictions of Gibbs’ canonical ensemble. This variety could be extended 
further by including one or more relaxation times. To choose among the combinations of 
{ a,/3 } it is convenient to use computerized searches, either seeking minimum deviations 
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Figure 4: ( a,/3 ) = ( 0.411,0.689 ) . A close inspection shows 26 prominent “holes” with 

a total measure less than one-half percent of the total. The Figure also illustrates the typical 
figure-eight-shaped nullcline where the trajectory motion is parallel to the C = 0 plane. Because 
the deviations of the various one-dimensional probability densities are visually indistinguishable 
from Gibbs’ canonical distributions none of them is shown in Figures 4-6. Here Ai = 0.1621 in 
the chaotic sea. 

with Monte Carlo searches^^ or by choosing minima from grid-based arrays of { a, P ) 
resnlts. 

We nse standard Rnnge-Kutta integration methods thronghont this work, fonrth-order 
and fifth-order as well as two types of “adaptive” integrators. In the adaptive cases the 
timestep dt is donbled if the error is “small” ( typically 10“^® ) and halved if the error 
is “large” ( typically 10“^^ ). The error estimate compares either fonrth-order and hfth- 
order resnlts with the same timestep dt or two fonrth-order resnlts with dt and dt/2 . 
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Figure 5: { a, (5 ) = { 0.354,0.746 ) likewise provides “reasonable” distributions and moments, 
but has four holes where nested tori penetrate the cross section. The tenfold “zoom” of one hole 
shows their roughly triangular shape. Here Ai = 0.1525 in the chaotic sea. 

There are no significant differences in the conclusions reached with any of these several 
methods. We carried out independent calculations in Nevada and in Wisconsin. 

With either two-parameter approach, the computation of moments is straightforward. 
Optimizing the dynamics so as to seek out the Gibbs distribution can be accomplished by 
evaluating either [ 1 ] moments of the distribution f{q,p) or [ 2 ] values of the distribution 
itself. To follow the hrst approach we minimized the summed-up squared deviations of 
the first five nonvanishing Gibbsian moments : 

^ [ (qVT^) - 3 ] + [ {qV/T^) - 1 ] + [ (pVT^) - 3 ] + 

[ (gVT) - 1 ] + [ (pVT) - 1 ] . 
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Figure 6: («,/?) = ( 0.273,0.827 ) is just one of an infinite number of combinations that is 

apparently ergodic. Ai = 0.1450 . 

We evaluated cr^ for thousands of runs. Each used 100 million timesteps for a particular 
pair of candidate values ( a,f3 ). Time-averaged values of suggest the range 0.25 < 
a < 0.65 with 1 < a + (3 < 1.1 . Figures 4-6 shows three typical cross-sectional plots of 
{ } sections selected in this way. Figure 4, with ( a,/5 ) = ( 0.411,0.689 ) , shows 

26 noticeable “holes” with a total measure near one percent of the total. The Figure also 
illustrates the hgure-eight-shaped nullcline where trajectories move parallel to the C = 0 
plane with ^ = 0 ; 

/?[ (gVT) -!]+«[ (pVT)2 - 3(pVT) ] ^ 0 . 

Despite all the holes indicating nested tori the distribution functions and their Erst several 
moments are close to the Gibbsian ergodic values. Evidently there is no real substitute 
for looking at the sections themselves. 

Figure 5, with ( q;,/9 ) = ( 0.354,0.746 ) , likewise provides “reasonable” distributions 
and moments, but has four holes where nested tori penetrate the cross section. Figures 
4 and 5 hint at the extensive zoo of topologies hidden in the { a, (3 ) plane. There are in 
addition patches of values which evidently correspond to ergodicity. Two examples which 
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we think are likely ergodic are : 


( 0.273,0.827 ) and ( 0.274,0.826 ) . 

The hrst of these is plotted in Figure 6 . It has no noticeable “holes”, indicating ergodicity 
within an accuracy of about one part per hundred thousand. Again the Figure-Eight- 
Shaped white space indicates the nullcline, which depends weakly on the precise value of 
the ratio {a/(3) . 

One can only search for “holes” in sections visually. An example, which we thought 
to be ergodic after a cursory inspection, is the combination ( a,/? ) = ( 0.495,0.555 ) , 
not shown here because visualizing the holes requires magnihcation. A close inspection 
of the C = 0 section reveals 36 tiny holes (!) corresponding to a single thin set of nested 
tori, including two near { q,p ) = { ±1.5, 0.0 ) . These tori cross the C = 0 plane in 36 
separate places. 

Because the numerical value of the summed squared-moment errors depends upon both 
the initial conditions and the length of the trajectories a reasonable procedure is to inves¬ 
tigate visually, as second and third criteria for ergodicity, the distributions themselves as 
well as their {q,p) cross sections. Such inspections reveal the {a, P ) pairs most promising 
from the standpoint of ergodicity. The distributions found for the cross sections of the 
Figures 1-3 are included in those Figures. From the visual standpoint such histograms 
show no signihcant deviations from the ergodic distribution in the data displayed for £g- 
ures 4 and 5, despite the clear nonergodicity seen in the cross sections. Our results suggest 
overall that a visual inspection of two-dimensional cross sections is the most reliable way 
to identify ergodicity in these three-dimensional dynamical systems. 

An alternative method for evaluating ergodicity is to compute the probability that a 
measured distribution of data points, such as { g* } or { pj } or { } comes from the 

expected Gaussian distribution of such points. Comparing the probabilities for the three 
choices demonstrates the accuracy of such a test. So long as the sampling bins contain 
many points the mean-squared deviation of the bin populations should be approximately 
equal to Gibbs’ value. Such tests implementing criteria can serve as useful indicators 
for deviations from ergodicity. In any doubtful case visual inspection is the only reliable 
criterion. 
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IV. PHASE-SPACE DENSITY FLOWS 


An apparent alternative to solving the motion equations { C } = for a specimen 
oscillator is the solution of Liouville’s continuity equation, /// = —V • n , so as to 
study the details of the convergence ( or lack of it ) to Gibbs’ canonical distribution. We 
briefly considered this approach and developed a straightforward hnite-difference program 
simulating the three-dimensional flow of the probability density f{q, p, t) . This program 
quickly led to negative densities. A conservative approach, passing probabilities between 
adjacent cells, can be implemented with a swarm of N moving particles, all of equal 
probability. The instantaneous summed-up density of these particles at any point in 
phase space can be made continuous and twice-differentiable by dehning and computing a 
smooth-particle density. This idea is simplest to implement using Lucy’s weight function^^ 
w{r < h) with a range h of order two or three times the nearest-neighbor particle spacing: 

N 

f{q,p, C) = '^2w{r — ri) ; w{z = r/h < 1) oc (1 + 32;)(1 — z)^ . 

i 

We explored this idea using an initial condition /(O < g,p, C < 1) = 1 and noticed 
that such a focalized initial value requires several Lyapunov times to smooth out. The 
particulate basis of the density guarantees that there is no tendency for this solution of the 
Liouville flow to stabilize. The time-reversible nature of the flow guarantees that a smooth 
stationary solution can only be obtained by adding in a time-averaging step. A detailed 
investigation of these ideas is likely worthwhile in that the compact three-dimensional 
nature of these flows makes visualization easy. 


V. SUMMARY 

It appears highly likely that a single thermostat variable is enough to provide Gibbs’ 
canonical distribution for a thermostated harmonic oscillator. This question has stimu¬ 
lated a relatively complex and varied literature over a 30-year period. The mathematicians 
are content to prove nonergodicity^. The computational physicists, ourselves included, 
have been prone to give up on the possibility of ergodicity with a single thermostat 
variable^. Thus our hnding that one thermostat is enough was a pleasant surprise. The 
{ a, P ) model detailed here seems likely to be the simplest smoothly deterministic, er- 
godic, time-reversible, and chaotic system for which the phase-space distribution is exactly 
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known. 
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